Preparations

library(caret) # load this first as it loads plyr
library(gbm)
library(car)
library(grid)
library(gridExtra)
library(randomForest)
library(xtable)
library(kableExtra)
library(abt)
library(tidyverse)
source('FISH_functions.R')

Raw fish data

fish = read.csv('data/Fish benthic physical sitelevel 2019.csv', strip.white=TRUE)
fish %>% glimpse
## Observations: 619
## Variables: 300
## $ YEAR                        <int> 2007, 2007, 2007, 2007, 2007, 2007, …
## $ REGION                      <fct> Keppel, Keppel, Keppel, Keppel, Kepp…
## $ EXPOSURE                    <fct> Semi-Exposed, Exposed, Semi-Exposed,…
## $ NTR                         <fct> Fished, Fished, Fished, Fished, Fish…
## $ NTR.Pooled                  <fct> Fished, Fished, Fished, Fished, Fish…
## $ SITE                        <fct> GK1, GK2, GK3, GK4, GK7, GK8, H3, MI…
## $ lut.arge                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.carp                    <dbl> 13.3320, 11.3322, 4.6662, 8.6658, 16…
## $ lut.flma                    <dbl> 0.0000, 3.3330, 0.0000, 0.0000, 7.33…
## $ lut.fulv                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.kas                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.lemn                    <dbl> 0.0000, 0.6666, 0.0000, 0.0000, 0.00…
## $ lut.lutj                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.mono                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.quin                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ lut.rivu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.russ                    <dbl> 67.3266, 4.6662, 0.0000, 5.3328, 116…
## $ lut.seba                    <dbl> 1.9998, 0.6666, 9.9990, 0.0000, 0.00…
## $ lut.vitt                    <dbl> 377.9622, 10.6656, 16.6650, 9.3324, …
## $ mac.macu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pagrus                      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sym.nema                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.barb                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ par.boides                  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.bifa                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.indi                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ par.multi                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.cili                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ upe.trag                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ platax                      <dbl> 0.0000, 0.6666, 0.0000, 0.6666, 0.66…
## $ Platycephalus.spp.          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.atki                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ let.hara                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.lati                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ let.lent                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ let.mini                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.nebu                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ let.oliv                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.obso                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.orna                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ gym.spp                     <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ mon.gran                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sco.bili                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sco.marg                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ sco.mono                    <dbl> 4.6662, 1.9998, 3.9996, 5.3328, 3.99…
## $ sco.line                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ any.leuc                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ath.roga                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cep.argu                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cep.boen                    <dbl> 21.9978, 5.9994, 7.3326, 8.6658, 3.3…
## $ cep.cyan                    <dbl> 1.3332, 0.0000, 1.9998, 1.3332, 0.00…
## $ cep.micr                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cro.alti                    <dbl> 0.6666, 0.0000, 0.6666, 0.6666, 0.00…
## $ epi.caer                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.coio                    <dbl> 0.0000, 0.0000, 0.0000, 1.3332, 0.00…
## $ epi.cora                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.fasc                    <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 0.00…
## $ epi.fusc                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.howl                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.lanc                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ epi.hexa                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.merr                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.ongu                    <dbl> 0.6666, 0.0000, 0.6666, 0.0000, 0.00…
## $ epi.quoy                    <dbl> 6.6660, 2.6664, 12.6654, 1.9998, 1.9…
## $ epi.sexf                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pms.laev                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pms.leop                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.66…
## $ pms.macu                    <dbl> 5.9994, 3.9996, 17.3316, 7.3326, 9.9…
## $ pte.ante                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pte.voli                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ dip.bifa                    <dbl> 1.9998, 1.9998, 7.9992, 1.9998, 1.99…
## $ dia.pict                    <dbl> 23.3310, 3.3330, 4.6662, 1.9998, 9.9…
## $ ple.chae                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ ple.flav                    <dbl> 1.9998, 1.3332, 2.6664, 0.6666, 0.00…
## $ ple.gibb                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.less                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.line                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.picu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.unic                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ bodianus                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ oxy.diag                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ che.chlo                    <dbl> 1.3332, 0.6666, 5.9994, 0.6666, 1.33…
## $ che.fasc                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ che.tril                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ che.undu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cho.anch                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cho.ceph                    <dbl> 0.6666, 0.0000, 0.6666, 1.3332, 0.66…
## $ cho.cya                     <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cho.fasc                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cho.grap                    <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 1.33…
## $ cho.scho                    <dbl> 1.3332, 0.0000, 0.6666, 0.6666, 0.00…
## $ cho.vitt                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ epb.insi                    <dbl> 0.6666, 0.0000, 0.0000, 0.0000, 1.33…
## $ gom.vari                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ hem.melt                    <dbl> 17.3316, 9.3324, 25.9974, 13.9986, 2…
## $ hem.fasc                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ psa.wai                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.bloc                    <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 4.66…
## $ aca.duss                    <dbl> 2.6664, 0.0000, 0.0000, 0.0000, 9.99…
## $ aca.gram                    <dbl> 0.0000, 0.0000, 0.6666, 1.3332, 4.66…
## $ aca.line                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.nans                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.ncus                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.nuda                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.xant                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ aca.thom                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.spp                     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cte.bino                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cte.stri                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.annu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.brach                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.brevi                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.litu                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.tong                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.unic                    <dbl> 0.0000, 1.9998, 0.0000, 0.0000, 0.00…
## $ pri.spp                     <dbl> 0.0000, 0.0000, 0.0000, 3.9996, 18.6…
## $ kyp.spp                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ zeb.scop                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ zeb.veli                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ zan.corn                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ bol.muri                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cal.caro                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cet.bico                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chl.blee                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chl.micr                    <dbl> 0.6666, 0.0000, 0.0000, 0.0000, 0.00…
## $ chs.sord                    <dbl> 10.6656, 2.6664, 0.0000, 0.0000, 0.0…
## $ hip.long                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.alti                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.cham                    <dbl> 6.6660, 0.6666, 0.0000, 0.0000, 0.00…
## $ sca.dimi                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.flav                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.fors                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.fren                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ sca.ghob                    <dbl> 7.9992, 3.9996, 1.9998, 0.6666, 13.3…
## $ sca.glob                    <dbl> 0.6666, 7.9992, 0.0000, 0.0000, 0.66…
## $ sca.nigr                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.ovic                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.psit                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.rivu                    <dbl> 37.3296, 25.3308, 19.3314, 7.9992, 4…
## $ sca.rubr                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.schl                    <dbl> 0.6666, 0.0000, 0.0000, 0.0000, 0.00…
## $ sca.spin                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.tric                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.spp                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.arge                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.cana                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.cora                    <dbl> 0.0000, 1.3332, 0.0000, 0.0000, 0.00…
## $ sig.doli                    <dbl> 5.3328, 3.9996, 1.3332, 2.6664, 6.66…
## $ sig.fusc                    <dbl> 0.0000, 0.0000, 11.9988, 6.6660, 0.0…
## $ sig.javu                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.line                    <dbl> 0.0000, 0.0000, 0.0000, 1.3332, 0.00…
## $ sig.puel                    <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 0.00…
## $ sig.ptus                    <dbl> 0.0000, 0.0000, 0.0000, 1.9998, 0.00…
## $ sig..ptiss                  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.stell                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.vulp                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ mon.arge                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.auri                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cha.afas                    <dbl> 63.3270, 89.9910, 39.3294, 57.3276, …
## $ cha.baro                    <dbl> 0.0000, 1.3332, 0.0000, 0.0000, 0.00…
## $ cha.benn                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.ephi                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cha.citr                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.flav                    <dbl> 0.0000, 0.0000, 0.0000, 0.6666, 0.00…
## $ cha.line                    <dbl> 5.9994, 5.3328, 0.0000, 1.3332, 5.33…
## $ cha.lunu                    <dbl> 1.3332, 0.0000, 0.6666, 1.3332, 0.00…
## $ cha.ttus                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.mela                    <dbl> 4.6662, 1.3332, 1.9998, 0.0000, 0.66…
## $ cha.ocel                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.orna                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.pleb                    <dbl> 0.0000, 0.0000, 0.6666, 0.6666, 0.00…
## $ cha.rain                    <dbl> 3.3330, 13.9986, 1.3332, 3.3330, 3.9…
## $ cha.raff                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.seme                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.spec                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cha.tlis                    <dbl> 1.3332, 2.6664, 0.0000, 0.0000, 0.00…
## $ cha.ulie                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cha.vaga                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chm.rost                    <dbl> 7.3326, 3.3330, 1.9998, 1.3332, 6.66…
## $ cor.alti                    <dbl> 3.9996, 1.9998, 0.0000, 0.6666, 1.99…
## $ cor.chry                    <dbl> 0.0000, 1.3332, 0.0000, 0.6666, 1.99…
## $ hen.acum                    <dbl> 0.0000, 1.3332, 0.0000, 0.0000, 0.00…
## $ hen.mono                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ hen.vari                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ mct.stri                    <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ prc.ocel                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.bico                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.bisp                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.nox                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.tibi                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.vrol                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chd.doub                    <dbl> 5.9994, 3.9996, 0.6666, 1.3332, 7.99…
## $ chd.mere                    <dbl> 5.3328, 5.3328, 1.9998, 2.6664, 4.66…
## $ pmc.impe                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pmc.semi                    <dbl> 0.0000, 3.3330, 0.0000, 0.0000, 0.66…
## $ pmc.sext                    <dbl> 1.3332, 1.3332, 0.0000, 0.6666, 0.66…
## $ pmc.xant                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pyg.diac                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ abu.spp                     <int> 406, 48, 0, 0, 238, 0, 30, 112, 38, …
## $ acn.poly                    <int> 90, 48, 74, 58, 124, 22, 152, 78, 62…
## $ amb.aure                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ amb.cura                    <int> 0, 0, 4, 2, 0, 0, 0, 0, 0, 2, 0, 0, …
## $ amb.leuc                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ amp.spp                     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ amp.peri                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pre.spp                     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chl.labi                    <int> 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.ambo                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.alis                    <int> 60, 0, 12, 0, 0, 0, 138, 34, 16, 64,…
## $ chr.apes                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.marg                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.niti                    <int> 8920, 9680, 7340, 15840, 12210, 240,…
## $ chr.retr                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.tern                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.webe                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.xant                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chy.roll                    <int> 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chy.talb                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chy.flav                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ das.spp                     <int> 12, 0, 0, 0, 0, 0, 0, 0, 0, 24, 4, 0…
## $ das.reti                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ das.trim                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ dis.spp                     <int> 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ hgy.plag                    <int> 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0,…
## $ neg.mela                    <int> 14, 6, 2, 4, 12, 0, 4, 0, 8, 0, 0, 2…
## $ parma                       <int> 10, 14, 0, 0, 4, 2, 2, 2, 4, 0, 4, 0…
## $ neg.nigr                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pgy.dick                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pgy.lacr                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.adel                    <int> 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, …
## $ pom.ambo                    <int> 0, 0, 8, 0, 0, 0, 18, 6, 4, 32, 0, 0…
## $ pom.aust                    <int> 58, 140, 14, 54, 28, 0, 216, 156, 21…
## $ pom.bank                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.brac                    <int> 84, 0, 4, 6, 52, 0, 90, 14, 10, 0, 6…
## $ pom.chry                    <int> 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, …
## $ pom.coel                    <int> 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.gram                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.lepi                    <int> 8, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.molu                    <int> 602, 294, 426, 634, 542, 190, 564, 1…
## $ pom.naga                    <int> 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.phil                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.reid                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.trip                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.vaiu                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.ward                    <int> 258, 224, 206, 338, 200, 456, 136, 1…
## $ ste.apic                    <int> 2, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, …
## $ ste.fasc                    <int> 8, 8, 4, 10, 8, 12, 0, 4, 8, 0, 4, 6…
## $ ste.livi                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ste.nigr                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ oxy.long                    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Anampses                    <int> 0, 6, 0, 2, 8, 0, 0, 0, 0, 0, 6, 2, …
## $ Coris                       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Halichoeres                 <int> 14, 10, 20, 30, 8, 8, 26, 6, 10, 48,…
## $ Labroides                   <int> 18, 10, 10, 12, 16, 4, 10, 10, 4, 2,…
## $ Labropsis                   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Psuedolabrus.guentheri      <int> 20, 58, 20, 24, 26, 104, 16, 22, 20,…
## $ Stethojoulis                <int> 12, 4, 16, 28, 2, 30, 16, 2, 0, 84, …
## $ Thalasoma                   <int> 24, 8, 12, 12, 24, 2, 30, 10, 16, 18…
## $ Macropharyngodon            <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Labrichthys                 <int> 8, 0, 0, 6, 0, 2, 0, 2, 4, 4, 4, 0, …
## $ Total.fish.density          <dbl> 11381.9272, 10805.9752, 8401.3122, 1…
## $ Total.fish.species.richness <dbl> 35.2, 29.4, 27.6, 26.4, 29.2, 17.8, …
## $ Prey.density                <int> 10654, 10558, 8190, 17060, 13502, 10…
## $ Prey.biomass                <dbl> 44.760041, 30.942239, 23.789262, 46.…
## $ Plectropomus.total.density  <dbl> 5.9994, 3.9996, 17.3316, 7.3326, 10.…
## $ Plectropomus.total.biomass  <dbl> 3.980542, 5.371889, 9.546584, 8.9263…
## $ Plectropomus.legal.density  <dbl> 1.9998, 2.6664, 4.6662, 3.9996, 1.99…
## $ Plectropomus.legal.biomass  <dbl> 2.4194697, 4.5827655, 5.0227904, 7.1…
## $ BE                          <dbl> 121.9966, 119.9976, 114.6630, 129.33…
## $ GRAZ                        <dbl> 72.6594, 47.9952, 36.6630, 26.6640, …
## $ COR                         <dbl> 103.3246, 119.3214, 45.9954, 71.9934…
## $ OM                          <dbl> 499.9996, 99.9996, 74.0000, 59.9998,…
## $ PL                          <dbl> 9000, 9680, 7366, 15842, 12210, 240,…
## $ CA                          <dbl> 521.9478, 45.9954, 65.9934, 43.9956,…
## $ PI                          <dbl> 7.9992, 6.6660, 25.9974, 10.6656, 12…
## $ FA                          <int> 290, 246, 218, 348, 212, 474, 138, 1…
## $ slope                       <dbl> 2.60, 2.64, 2.00, 2.40, 1.96, 2.20, …
## $ rugosity                    <dbl> 2.68, 3.08, 3.00, 3.00, 3.00, 2.68, …
## $ SCI                         <dbl> 7.040, 8.088, 6.000, 7.200, 5.880, 5…
## $ LCC_.                       <dbl> 60.0, 78.8, 77.6, 77.6, 57.2, 13.6, …
## $ LHC_.                       <dbl> 56.4, 68.8, 77.6, 75.6, 55.6, 13.6, …
## $ SC_.                        <dbl> 3.6, 10.0, 0.0, 2.0, 1.6, 0.0, 0.0, …
## $ MAC_.                       <dbl> 2.4, 0.4, 9.6, 4.4, 3.6, 4.0, 29.2, …
## $ SPO_.                       <dbl> 1.6, 0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0…
## $ Turf_.                      <dbl> 26.0, 17.6, 4.4, 11.6, 28.0, 71.2, 6…
## $ Unconsolidated_.            <dbl> 9.6, 3.2, 8.4, 6.0, 11.2, 11.2, 7.2,…
## $ Benthic.richness            <dbl> 6.6, 5.8, 2.2, 6.0, 3.8, 2.0, 3.4, 5…
## $ Coral_Morph.Diversity       <dbl> 5.4, 5.6, 1.2, 5.0, 3.4, 1.4, 2.4, 4…
## $ ChlA                        <dbl> 0.45221510, 0.45221510, 0.41940367, …
## $ kd490                       <dbl> 0.06519040, 0.06519040, 0.06261333, …
## $ SSTmean                     <fct> 23.80386083, 23.80357108, 23.7912496…
## $ SSTanom                     <dbl> -0.3130472, -0.3140394, -0.3159858, …
## $ wave.exposure.index         <dbl> 0.007256172, 0.160425811, 0.00099007…
## $ Corrected.depth             <dbl> 4.03, 6.39, 5.12, 5.47, 3.73, 3.02, …
## $ maxDHW                      <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
## $ deltaLHC                    <dbl> 14.0, 7.6, 6.4, 16.8, 8.0, 2.8, 4.8,…
## $ deltaMA                     <dbl> -6.8, -1.2, -3.6, 0.4, -6.4, -41.2, …
## $ Cyclone                     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ LT.Fprimary                 <dbl> 0.087, 0.087, 0.087, 0.031, 0.025, 0…
## $ Exposure.to.primary.weeks   <int> 2, 2, 2, 2, 3, 2, 3, 1, 1, 0, 3, 2, …
var.lookup = rbind(
    data.frame(pretty.name='Total density', Field.name='Total.fish.density', Abbreviation='TFD', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Total species richness', Field.name='Total.fish.species.richness', Abbreviation='TFSR', Family='gaussian', Type='Response', Transform='I', Groupby=''),
    data.frame(pretty.name='Benthic invertivores', Field.name='BE', Abbreviation='BE', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Grazers', Field.name='GRAZ', Abbreviation='GRAZ', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Corallivores', Field.name='COR', Abbreviation='COR', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Omnivores', Field.name='OM', Abbreviation='OM', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Planktivores', Field.name='PL', Abbreviation='PL', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Carnivores', Field.name='CA', Abbreviation='CA', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Piscivores', Field.name='PI', Abbreviation='PI', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Farmers', Field.name='FA', Abbreviation='FA', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Plectropomus total density', Field.name='Plectropomus.total.density', Abbreviation='PTD', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Plectropomus total biomass', Field.name='Plectropomus.total.biomass', Abbreviation='PTB', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Plectropomus legal density', Field.name='Plectropomus.legal.density', Abbreviation='PLD', Family='gaussian', Type='Response', Transform='log', Groupby=''),
    data.frame(pretty.name='Plectropomus legal biomass', Field.name='Plectropomus.legal.biomass', Abbreviation='PLB', Family='gaussian', Type='Response', Transform='log', Groupby=''),

    data.frame(pretty.name='Region', Field.name='REGION', Abbreviation='REGION', Family=NA, Type='Predictor', Transform='I', Groupby=''),
    data.frame(pretty.name='NTR Pooled', Field.name='NTR.Pooled', Abbreviation='NTR.Pooled', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='LHC % (live hard coral cover)', Field.name='LHC_.', Abbreviation='LHC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='SC % (soft coral)', Field.name='SC_.', Abbreviation='SC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='MA % (macroalgal cover)', Field.name='MAC_.', Abbreviation='MA', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Turf %', Field.name='Turf_.', Abbreviation='TURF', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Unconsolidated %', Field.name='Unconsolidated_.', Abbreviation='UC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Benthic richness', Field.name='Benthic.richness', Abbreviation='BR', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Coral morphological diversity', Field.name='Coral_Morph.Diversity', Abbreviation='CMD', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Slope', Field.name='slope', Abbreviation='SLOPE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Rugosity', Field.name='rugosity', Abbreviation='RUG', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='SCI (structural complexity)', Field.name='SCI', Abbreviation='SCI', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Chla', Field.name='ChlA', Abbreviation='CHL', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Kd490', Field.name='kd490', Abbreviation='KD490', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='SST mean', Field.name='SSTmean', Abbreviation='SSTMEAN', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='SST anom', Field.name='SSTanom', Abbreviation='SSTANOM', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Wave exposure', Field.name='wave.exposure.index', Abbreviation='WAVE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Corrected depth', Field.name='Corrected.depth', Abbreviation='DEPTH', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Max DHW', Field.name='maxDHW', Abbreviation='DHW', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Cyclone', Field.name='Cyclone', Abbreviation='CYCLONE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Exposure to primary weeks', Field.name='Exposure.to.primary.weeks', Abbreviation='EXP', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Prey density', Field.name='Prey.density', Abbreviation='PREY.DENSITY', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
    data.frame(pretty.name='Prey biomass', Field.name='Prey.biomass', Abbreviation='PREY.BIOMASS', Family=NA, Type='Predictor', Transform='I', Groupby='Region')
)
save(var.lookup, file='data/var.lookup.RData') 
names = with(var.lookup, setNames(as.character(Field.name), Abbreviation))
## exclude those whose names equal their values otherwise there will be duplicate fields created in the fish data
names = names[names(names)!=names]
## Create duplicates of the fields that are not already abbreviated
fish = fish %>%
    mutate(SSTmean=ifelse(SSTmean=='#N/A',NA,SSTmean)) %>%
    mutate_at(as.character(var.lookup$Field.name), list(A=~I)) %>%
    rename(!!! gsub('(.*)','\\1_A',names)) %>%
    mutate(REGION=factor(REGION, levels=c('Palm','Magnetic','Whitsunday','Keppel')))
save(fish, file='data/fish.RData')

Response Variables:

Variable Data field Abbreviation Distribution
Total density Total.fish.density TFD Gaussian (log)
Total species richness Total.fish.species.richness TFSR Gaussian
Benthic invertivores BE BE Gaussian (log)
Grazers GRAZ GRAZ Gaussian (log)
Corallivores COR COR Gaussian (log)
Omnivores OM OM Gaussian (log)
Planktivores PL PL Gaussian (log)
Carnivores CA CA Gaussian (log)
Piscivores PI PI Gaussian (log)
Farmers FA FA Gaussian (log)
Plectropomus total density Plectropomus.total.density PTD Gaussian (log)
Plectropomus total biomass Plectropomus.total.biomass PTB Gaussian (log)
Plectropomus legal density Plectropomus.legal.density PLD Gaussian (log)
Plectropomus legal biomass Plectropomus.legal.biomass PLB Gaussian (log)

Exploratory data analyses

EDA_histograms = function(var='', group='',dat=NULL, var.lookup) {
    pretty.name = as.character(var.lookup$pretty.name[var.lookup$Abbreviation==var])
    if (is.numeric(dat[,var])) {
        grid.arrange(
            ggplot(dat) + geom_histogram(aes(x=!!sym(var))) +
            scale_x_continuous(pretty.name),
            ggplot(dat) + geom_histogram(aes(x=!!sym(var))) +
            scale_x_log10(pretty.name),
            nrow=1
        )
    } else {
        grid.arrange(
            ggplot(dat) + geom_bar(aes(x=!!sym(var), fill=!!sym(group))) +
            scale_x_discrete(pretty.name),
            nrow=1
        )
    }
}
EDA_density = function(var='', group='', dat=NULL, var.lookup) {
    pretty.name = as.character(var.lookup$pretty.name[var.lookup$Abbreviation==var])
    if (is.numeric(dat[,var])) {
        grid.arrange(
            ggplot(dat) + geom_density(aes(x=!!sym(var), fill=!!sym(group)), alpha=0.5) +
            scale_x_continuous(pretty.name),
            ggplot(dat) + geom_density(aes(x=!!sym(var), fill=!!sym(group)), alpha=0.5) +
            scale_x_log10(pretty.name),
            nrow=1
        )
    } else {
    }
}

Total density

EDA_histograms(var='TFD', dat=fish, var.lookup=var.lookup)

EDA_density(var='TFD', group='REGION', dat=fish, var.lookup=var.lookup)

Total species richness

EDA_histograms(var='TFSR', dat=fish, var.lookup=var.lookup)

EDA_density(var='TFSR', group='REGION', dat=fish, var.lookup=var.lookup)

Benthic invertivores

EDA_histograms(var='BE', dat=fish, var.lookup=var.lookup)

EDA_density(var='BE', group='REGION', dat=fish, var.lookup=var.lookup)

Grazers

EDA_histograms(var='GRAZ', dat=fish, var.lookup=var.lookup)

EDA_density(var='GRAZ', group='REGION', dat=fish, var.lookup=var.lookup)

Corallivores

EDA_histograms(var='COR', dat=fish, var.lookup=var.lookup)

EDA_density(var='COR', group='REGION', dat=fish, var.lookup=var.lookup)

Omnivores

EDA_histograms(var='OM', dat=fish, var.lookup=var.lookup)

EDA_density(var='OM', group='REGION', dat=fish, var.lookup=var.lookup)

Planktivores

EDA_histograms(var='PL', dat=fish, var.lookup=var.lookup)

EDA_density(var='PL', group='REGION', dat=fish, var.lookup=var.lookup)

Carnivores

EDA_histograms(var='CA', dat=fish, var.lookup=var.lookup)

EDA_density(var='CA', group='REGION', dat=fish, var.lookup=var.lookup)

Piscivores

EDA_histograms(var='PI', dat=fish, var.lookup=var.lookup)

EDA_density(var='PI', group='REGION', dat=fish, var.lookup=var.lookup)

Farmers

EDA_histograms(var='FA', dat=fish, var.lookup=var.lookup)

EDA_density(var='FA', group='REGION', dat=fish, var.lookup=var.lookup)

Plectropomus total density

EDA_histograms(var='PTD', dat=fish, var.lookup=var.lookup)

EDA_density(var='PTD', group='REGION', dat=fish, var.lookup=var.lookup)

Plectropomus total biomass

EDA_histograms(var='PTB', dat=fish, var.lookup=var.lookup)

EDA_density(var='PTB', group='REGION', dat=fish, var.lookup=var.lookup)

Predictor Variables:

Variable Data field Abbreviation
Region REGION REGION
NTR Pooled NTR.Pooled NTR.Pooled
LHC % (live hard coral cover) LHC_. LHC
SC % (soft coral) SC_. SC
MA % (macroalgal cover) MAC_. MA
Turf % Turf_. TURF
Unconsolidated % Unconsolidated_. UC
Benthic richness Benthic.richness BR
Coral morphological diversity Coral_Morph.Diversity CMD
Slope slope SLOPE
Rugosity rugosity RUG
SCI (structural complexity) SCI SCI
Chla ChlA CHL
Kd490 kd490 KD490
SST mean SSTmean SSTMEAN
SST anom SSTanom SSTANOM
Wave exposure wave.exposure.index WAVE
Corrected depth Corrected.depth DEPTH
Max DHW maxDHW DHW
Cyclone Cyclone CYCLONE
Exposure to primary weeks Exposure.to.primary.weeks EXP

FOR Plectropomus response variables ONLY, add:

Variable Data field Abbreviation
Prey density (for Plectropomus density and legal density) Prey.density PREY.DENSITY
Prey biomass (for Plectropomus biomass and legal biomass) Prey.biomass PREY.BIOMASS

Species-specific analyses: list of species in the attached file. There are two lists there: the single-column list is for the initial analysis between regions, and the table is for the species lists to use for each region individually.

Raw selections data

selections = read.csv('data/Species selection.csv', strip.white=TRUE)
selections %>% glimpse
## Observations: 264
## Variables: 50
## $ Species                 <fct> aca.bloc, aca.duss, aca.gram, aca.line, …
## $ Palms                   <dbl> 2517.48, 3359.97, 1045.62, 0.00, 13.32, …
## $ Magnetic                <dbl> 579.42, 679.32, 126.54, 0.00, 0.00, 0.00…
## $ Whitsundays             <dbl> 1728.27, 965.70, 3716.28, 13.32, 76.59, …
## $ Keppels                 <dbl> 789.21, 346.32, 1312.02, 0.00, 23.31, 0.…
## $ Palms...total           <dbl> 0.150705201, 0.201139613, 0.062594488, 0…
## $ Magnetic...total        <dbl> 0.94095997, 1.10319444, 0.20549700, 0.00…
## $ Whitsundays...total     <dbl> 0.053913389, 0.030125015, 0.115929368, 0…
## $ Keppels...total         <dbl> 0.0145102620, 0.0063673720, 0.0241225460…
## $ Palms...FG              <dbl> 2.816481633, 3.759034349, 1.169808509, 0…
## $ Magnetic...FG           <dbl> 12.3404255, 14.4680851, 2.6950355, 0.000…
## $ Whitsundays...FG        <dbl> 1.364317447, 0.762335375, 2.933676822, 0…
## $ Keppels...FG            <dbl> 1.697099893, 0.744718940, 2.821339062, 0…
## $ X                       <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ X.1                     <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ Species.1               <fct> pom.molu, pom.brac, pom.ward, acn.poly, …
## $ Palms...total.1         <dbl> 23.9783308, 16.8312258, 7.5841087, 5.090…
## $ Species.2               <fct> pom.ward, pom.molu, Halichoeres, lut.flm…
## $ Magnetic...total.1      <dbl> 38.3094226, 11.1891446, 7.4864959, 6.451…
## $ Species.3               <fct> pom.molu, pom.brac, chr.niti, acn.poly, …
## $ Whitsundays...total.1   <dbl> 21.6939414, 15.6405509, 15.1548448, 8.69…
## $ Species.4               <fct> chr.niti, pom.ward, pom.molu, pom.aust, …
## $ Keppels...total.1       <dbl> 85.50374438, 4.50911903, 2.86450864, 1.4…
## $ X.2                     <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ X.3                     <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ Species.5               <fct> sca.rivu, sig.doli, sca.flav, chs.sord, …
## $ Palms...FG.1            <dbl> 37.6387750, 19.0187020, 7.9837568, 4.157…
## $ Species.6               <fct> sca.rivu, sca.ghob, aca.duss, aca.bloc, …
## $ Magnetic...FG.1         <dbl> 39.92907801, 17.23404255, 14.46808511, 1…
## $ Species.7               <fct> sca.rivu, sig.doli, sca.flav, sca.nigr, …
## $ Whitsundays...FG.1      <dbl> 31.6789779, 13.0411924, 12.1106175, 10.3…
## $ Species.8               <fct> sca.rivu, sca.ghob, sig.arge, sig.fusc, …
## $ Keppels...FG.1          <dbl> 47.26817043, 15.33834586, 11.57178661, 6…
## $ X.4                     <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ X.5                     <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ X.6                     <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ X5.most.abundant        <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ FG                      <fct> Grazers, , , , , Planktivores, , , , , C…
## $ Palms.1                 <fct> sca.rivu, sig.doli, sca.flav, chs.sord, …
## $ Magnetic.1              <fct> sca.rivu, sca.ghob, aca.duss, aca.bloc, …
## $ Whitsundays.1           <fct> sca.rivu, sig.doli, sca.flav, sca.nigr, …
## $ Keppels.1               <fct> sca.rivu, sca.ghob, sig.arge, sig.fusc, …
## $ X.7                     <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ X.8                     <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ X75..cumulative.minimum <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ FG.1                    <fct> Grazers, , , , , , , Planktivores, , , ,…
## $ Palms.2                 <fct> sca.rivu, sig.doli, sca.flav, chs.sord, …
## $ Magnetic.2              <fct> sca.rivu, sca.ghob, aca.duss, aca.bloc, …
## $ Whitsundays.2           <fct> sca.rivu, sig.doli, sca.flav, sca.nigr, …
## $ Keppels.2               <fct> sca.rivu, sca.ghob, sig.arge, sig.fusc, …

Exploratory data analyses

Total density

EDA_histograms(var='REGION', dat=fish, var.lookup=var.lookup)

NTR Polled

EDA_histograms(var='NTR.Pooled', dat=fish, var.lookup=var.lookup)

EDA_histograms(var='NTR.Pooled', group='REGION',dat=fish, var.lookup=var.lookup)

LHC % (live hard coral cover)

EDA_histograms(var='LHC', dat=fish, var.lookup=var.lookup)

EDA_density(var='LHC', group='REGION', dat=fish, var.lookup=var.lookup)

SC % (soft coral)

EDA_histograms(var='SC', dat=fish, var.lookup=var.lookup)

EDA_density(var='SC', group='REGION', dat=fish, var.lookup=var.lookup)

MA % (macroalgal cover)

EDA_histograms(var='MA', dat=fish, var.lookup=var.lookup)

EDA_density(var='MA', group='REGION', dat=fish, var.lookup=var.lookup)

Turf %

EDA_histograms(var='TURF', dat=fish, var.lookup=var.lookup)

EDA_density(var='TURF', group='REGION', dat=fish, var.lookup=var.lookup)

Unconsolidated %

EDA_histograms(var='UC', dat=fish, var.lookup=var.lookup)

EDA_density(var='UC', group='REGION', dat=fish, var.lookup=var.lookup)

Benthic richness

EDA_histograms(var='BR', dat=fish, var.lookup=var.lookup)

EDA_density(var='BR', group='REGION', dat=fish, var.lookup=var.lookup)

Coral morphological diversity

EDA_histograms(var='CMD', dat=fish, var.lookup=var.lookup)

EDA_density(var='CMD', group='REGION', dat=fish, var.lookup=var.lookup)

Slope

EDA_histograms(var='SLOPE', dat=fish, var.lookup=var.lookup)

EDA_density(var='SLOPE', group='REGION', dat=fish, var.lookup=var.lookup)

Rugosity

EDA_histograms(var='RUG', dat=fish, var.lookup=var.lookup)

EDA_density(var='RUG', group='REGION', dat=fish, var.lookup=var.lookup)

SCI (structural complexity)

EDA_histograms(var='SCI', dat=fish, var.lookup=var.lookup)

EDA_density(var='SCI', group='REGION', dat=fish, var.lookup=var.lookup)

Chlorophyl-a

EDA_histograms(var='CHL', dat=fish, var.lookup=var.lookup)

EDA_density(var='CHL', group='REGION', dat=fish, var.lookup=var.lookup)

Kd490

EDA_histograms(var='KD490', dat=fish, var.lookup=var.lookup)

EDA_density(var='KD490', group='REGION', dat=fish, var.lookup=var.lookup)

SST mean

EDA_histograms(var='SSTMEAN', dat=fish, var.lookup=var.lookup)

EDA_density(var='SSTMEAN', group='REGION', dat=fish, var.lookup=var.lookup)

SST anomoly

EDA_histograms(var='SSTANOM', dat=fish, var.lookup=var.lookup)

EDA_density(var='SSTANOM', group='REGION', dat=fish, var.lookup=var.lookup)

Wave exposure

EDA_histograms(var='WAVE', dat=fish, var.lookup=var.lookup)

EDA_density(var='WAVE', group='REGION', dat=fish, var.lookup=var.lookup)

Corrected Depth

EDA_histograms(var='DEPTH', dat=fish, var.lookup=var.lookup)

EDA_density(var='DEPTH', group='REGION', dat=fish, var.lookup=var.lookup)

Max DHW

EDA_histograms(var='DHW', dat=fish, var.lookup=var.lookup)

EDA_density(var='DHW', group='REGION', dat=fish, var.lookup=var.lookup)

Cyclone

EDA_histograms(var='CYCLONE', dat=fish, var.lookup=var.lookup)

EDA_density(var='CYCLONE', group='REGION', dat=fish, var.lookup=var.lookup)

Explore to primary water (weeks)

EDA_histograms(var='EXP', dat=fish, var.lookup=var.lookup)

EDA_density(var='EXP', group='REGION', dat=fish, var.lookup=var.lookup)

Results

formulas = list(
    all = ~REGION + NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
    Palm = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
    Magnetic = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
    Whitsunday = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
    Keppel = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP
)
load('data/fish.RData')
load('data/var.lookup.RData')
resp.lookup = var.lookup %>% filter(Type=='Response') %>% droplevels
pred.lookup = var.lookup %>% filter(Type=='Predictor') %>% droplevels

groupings.all = vector('list', nrow(pred.lookup))
names(groupings.all) = pred.lookup$Abbreviation
groupings.region = vector('list', nrow(pred.lookup))
names(groupings.region) = pred.lookup$Abbreviation
for (i in 1:nrow(pred.lookup)) {
    pred = as.character(pred.lookup[i,'Abbreviation'])
    groupings.all[[pred]] = ifelse(pred=='REGION',NA,'REGION')
    groupings.region[[pred]] = ifelse(pred=='NTR.Pooled',NA,'NTR.Pooled')
}
groupings.all = do.call('c',groupings.all)
groupings.region = do.call('c',groupings.region)

gfun = function(f) {
    form = attr(terms(formulas[[f]]),'term.labels')
    if (f=='all') return(groupings.all[form])
    else return(groupings.region[form])
}

analyses = vector('list',nrow(resp.lookup))
names(analyses) = resp.lookup$Abbreviation
for (i in 1:nrow(resp.lookup)) {
    resp = as.character(resp.lookup[i,'Abbreviation'])
    fun = as.character(resp.lookup[i,'Transform'])
    fam = as.character(resp.lookup[i,'Family'])
    
    analyses[[resp]] = list('formulas' = lapply(formulas, function(f) update(f, paste0(fun,'(',resp,') ~.'))),
                            'family' = fam)
    analyses[[resp]][['groups']] = sapply(c('all','Palm','Magnetic','Whitsunday','Keppel'), gfun, USE.NAMES = TRUE,simplify=FALSE)
}
for (a in 1:length(analyses)) {
    resp=names(analyses)[a]
    print(paste('Response =',resp))
    for (f in 1:length(analyses[[a]]$formulas)) {
        mod.name = names(analyses[[a]]$formulas)[f]
        print(paste('Model =',mod.name))
        MONOTONE = assignMonotone(fish, analyses[[a]]$formulas[[f]])
        if (mod.name=='all') {fish.sub=fish
        } else {fish.sub = fish %>% filter(REGION==mod.name)}
        if (any(fish.sub[,as.character(get_response(analyses[[a]]$formulas[[f]]))]==0)) {
            val = fish.sub[, as.character(get_response(analyses[[a]]$formulas[[f]]))]
            val=min(val[val>0], na.rm=TRUE)
            fish.sub[, as.character(get_response(analyses[[a]]$formulas[[f]]))] = fish.sub[,as.character(get_response(analyses[[a]]$formulas[[f]]))] + val
        }
        set.seed(123)
        #if ((analyses[[a]]$formulas[[f]])[[2]][[1]]=='log' & attr(terms(analyses[[a]]$formulas[[f]]), 'response')
        mod = abt(analyses[[a]]$formulas[[f]], data=fish.sub, distribution=analyses[[a]]$family,
                  cv.folds=10,interaction.depth=10,n.trees=10000, shrinkage=0.001, n.minobsinnode=2,
                  var.monotone=as.vector(MONOTONE))
        p=plot.abts(mod, var.lookup, center=FALSE, type='response', return.grid=TRUE,
                    groupby=na.omit(unique(analyses[[a]]$groups[[f]])))#
        thresholds=p$thresholds
        save(thresholds, file=paste0('data/thresholds_',mod.name,'_',resp,'.RData'))
        ps = p[['ps']]
        p = p[['p']]
        if (length(levels(fish$REGION))>1 || length(levels(fish$NTR.Pooled))>1) {
            p = common_legend(p)
            ps = common_legend(ps)
        }
        ## version for the supplimentary
        #do.call('grid.arrange', p) ## version for the supplimentary
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT.png'), do.call('grid.arrange', p), width=15, height=10, dpi=300)
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT.pdf'), do.call('grid.arrange', p), width=15, height=10)

        pred.1=stats.abt(mod, fitMethod=1, analysis = analyses[[a]]$groups[[f]])
        
        save(pred.1, file=paste0('data/pred.1_',mod.name,'_',resp,'.RData'))
        
        optim=summarize_values(pred.1$optim, type='optim') %>%
            dplyr::rename_at(vars(-contains('Var')), paste0, '.optim')
        
        rel.inf = summarize_values(pred.1$rel.imp, type='Rel.inf') %>%
            dplyr::rename_at(vars(-contains('Var')), paste0, '.rel.inf')
        
        R2=summarize_values(pred.1$R2.value, 'R2') %>%
            dplyr::rename_at(vars(-contains('Var')), paste0, '.R2')
        
        stats = optim %>% full_join(R2) %>% full_join(rel.inf) %>%
            left_join(var.lookup %>% dplyr::select(Var=Abbreviation, pretty.name))
        
        save(stats, file=paste0('data/stats_',mod.name,'_',resp,'.RData'))
        write.csv(stats %>% as.data.frame, file=paste0('output/data/stats_',mod.name,'_',resp,'.csv'), quote=FALSE, row.names=FALSE)

        ## Version full model with just the relative importance and the substantial panels
        ps1 = arrangeGrob(ps[[1]],ps[[length(ps)]], widths=c(3,1))
        numberOfPlots = length(ps)-2 #minus 1 since one of the items is the legend and minus one for the relative importance plot
        numberOfPlotRows = ceiling(numberOfPlots / 2)
        g=arrangeGrob(ps1, do.call('arrangeGrob', c(ps[c(-1, -length(ps))], list(ncol=2))), heights=c(2,numberOfPlotRows))
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT_short.png'), g, width=10, height=2+(1.7*numberOfPlotRows), dpi=300)
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT_short.pdf'), g, width=10, height=2+(1.7*numberOfPlotRows))

        ## Version with a grid of partials in the lower right corner of the influence figure
        ## The location and size of the figure in figure will be determined by the scale of relative influence
        ## ideally, we want the subfigure to be 3/4 of the width and 3/4 of the height
        ymax=max(pretty(as.numeric(as.character(stats$upper.rel.inf))))
        ymin=ymax*1/4
        xmax=(length(p)-2)*3/4
        if (ymin <= 100/(length(p)-2)) {  # if the left side of the subfigure is too close to the dashed vertical line
            ymax = ymax + 5
            ymin=ymax*1/4
            ps[[1]] = ps[[1]] + scale_y_continuous('Relative Importance', limits=c(0,ymax))
        }
        ps[2:(length(ps)-1)]=lapply(ps[2:(length(ps)-1)], function(f) f+theme(axis.title.y=element_blank())) 
        ps2=do.call('arrangeGrob', c(ps[c(length(ps),2:(length(ps)-1))], list(ncol=2)))
        #g= ps[[1]] + annotation_custom(grob=ps2, xmin=1, xmax=15, ymin=10, ymax=Inf)
        g= ps[[1]] + annotation_custom(grob=ps2, xmin=1, xmax=xmax, ymin=ymin, ymax=Inf)
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT_short_in.png'), g, width=10, height=8, dpi=300)
        ggsave(filename=paste0('output/figures/data.all.abt.',mod.name,'_',resp,'_ABT_short_in.pdf'), g, width=10, height=8)
                
        ## Version for vertical table of regional relative importance and substantial panels
        g=do.call('arrangeGrob', c(ps[1:4], list(ncol=1, heights=c(1.5, rep(1,3)))))
                                        #save(g, file=paste0('data/g.abt_',mod_num,'_',resp,'.RData'))
        
                                        #save(data.all.abt, file=paste0('data/data.all.abt_',mod_num,'_',resp,'.RData'))
        
                
        #save(rel.importance, file=paste0('data/rel.importance_',mod_num,'_',resp,'.RData'))
        #save(thresholds, file=paste0('data/thresholds_', mod_num,'_',resp,'.RData'))
                                        #write.table(bind_rows(thresholds, .id='Var'), file=paste0('data/thresholds_',mod_num,'_',resp,'.csv'), quote=FALSE, row.names=FALSE)
        rm(mod,pred.1,p,ps,ps1,g,stats,optim,rel.inf,R2,fit)
        gc()
    }
}


#summary(mod)
#plot(mod,c(15,1))
load('data/var.lookup.RData')
resp.lookup = var.lookup %>% filter(Type=='Response') %>% droplevels
mods = names(formulas)
for (a in 1:length(analyses)) {
    resp=names(analyses)[a]
    cat(paste('## ',resp.lookup$pretty.name[resp.lookup$Abbreviation==resp],' {.tabset .tabset-pills} \n\n'))
    for (f in 1:length(analyses[[a]]$formulas)) {
        mod.name = names(analyses[[a]]$formulas)[f]
        cat(paste('### ',mod.name,'\n\n'))
        cat(paste0('![](output/figures/data.all.abt.',mod.name,'_',resp,'_ABT_short_in.png)\n\n'))
        ## Now for a tables
        load(file=paste0('data/stats_',mod.name,'_',resp,'.RData'))
        pretty.stats(stats, var.lookup) %>% rmarkdown:::paged_table_html() %>% cat
        load(file=paste0('data/thresholds_',mod.name,'_',resp,'.RData'))
        pretty.thresholds(thresholds, stats, var.lookup) %>% rmarkdown:::paged_table_html() %>% cat
    }
}

Total density

all

Palm

Magnetic

Whitsunday

Keppel

Total species richness

all

Palm

Magnetic

Whitsunday

Keppel

Benthic invertivores

all

Palm

Magnetic

Whitsunday

Keppel

Grazers

all

Palm

Magnetic

Whitsunday

Keppel

Corallivores

all

Palm

Magnetic

Whitsunday

Keppel

Omnivores

all

Palm

Magnetic

Whitsunday

Keppel

Planktivores

all

Palm

Magnetic

Whitsunday

Keppel

Carnivores

all

Palm

Magnetic

Whitsunday

Keppel

Piscivores

all

Palm

Magnetic

Whitsunday

Keppel

Farmers

all

Palm

Magnetic

Whitsunday

Keppel

Plectropomus total density

all

Palm

Magnetic

Whitsunday

Keppel

Plectropomus total biomass

all

Palm

Magnetic

Whitsunday

Keppel